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ABSTRACT 

We study dynamics of bars in models of disk galaxies embeded in realistic dark matter 
halos. We find that disk thickness plays an important, if not dominant, role in the 
evolution and structure of the bars. We also make extensive numerical tests of different 
./V-body codes used to study bar dynamics. Models with thick disks typically used in 
this type of modeling (height-to-length ratio h z /R d = 0.2) produce slowly rotating, 
and very long, bars. In contrast, more realistic thin disks with the same parameters 
as in our Galaxy (h z /Rd ~ 0.1) produce bars with normal length i?bar ~ Rd, which 
rotate quickly with the ratio of the corotation radius to the bar radius TZ — 1.2 — 
1.4 compatible with observations. Bars in these models do not show a tendency to 
slow down, and may lose as little as 2-3 percent of their angular momentum due to 
dynamical friction with the dark matter over cosmological time. We attribute the 
differences between the models to a combined effect of high phase-space density and 
smaller Jeans mass in the thin disk models, which result in the formation of a dense 
central bulge. Special attention is paid to numerical effects such as the accuracy of 
orbital integration, force and mass resolution. Using three iV-body codes - Gadget, 
ART, and Pkdgrav - we find that numerical effects are very important and, if not 
carefully treated, may produce incorrect and misleading results. Once the simulations 
are performed with sufficiently small time-steps and with adequate force and mass 
resolution, all the codes produce nearly the same results: we do not find any systematic 
deviations between the results obtained with TREE codes (Gadget and Pkdgrav) and 
with the Adaptive-Mesh-Refinement (ART) code. 

Key words: galaxies:kinematics and dynamics — galaxies: evolution — galaxies: halos 
— methods:A— body simulations 



1 INTRODUCTION 



Barred galaxies re present a large fract i on (~65% ) of a ll 
spiral galaxies (e.g.. lEskridge et al.ll2000l ; ISheth et alj|2008h . 
Bars are ubiquitous. They are found in all type s of spirals: 
in large lenticular galaxies jAguerri et al.l|2005l), in normal 
spirals such as our Galaxy ([ Frcudcnreich 1998) and M31 



Beaton et al. 12007). and in 



|Athanassoula fc Beaton] 12006 
dwarf magellanic-type galaxies ( Valenzuela et all 20071 ) . An 
isolated stellar disk embeded into a dark matter halo spon- 
taneously forms a stellar bar as a r esult of the develop- 
ment of global disk instabilities (e.g., iBinnev fc Tremaine I 
1 19871 . Sec. 6.5). Bars continue to be closely scrutinized be- 
cause of their connection with the dark matter halo (e.g 



O'Neill fc Dubinski 1 120031: iHollev-Bockelmann et al. II2005I: 



Colm et al. 1120061 : lAthanassoula I I2007D . Because the bars 
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rotate inside massive dark matter halos, they lose some frac- 
tion of the angular momentum to their halos and tend to 
slow down with time l|Tremaine fc Weinberdll984l : IWeinberg| 
Il985l ). 

The formation of bars and associated pseudobulges is 
often consi dered as an alternative to the h ierarchical cluster- 
ing model (|Kormendv fc Kennicuttll2004l ). This appears to 
be incorrect: recent cosmological simulations indicate that 
the secular bulge formation is a part (not an alternative) 
of the hierarchical scenario. The simulations of the forma- 
tion of galaxies in the framework of the standard hierar- 
chical cosmological model indicate that bars form routinely 
in the course of as s embly of halos and galaxies inside them 
l|Maver et a l. 2008; Ceverinol l2008) . The simulations have a 
fine resolution of ~ 100 pc and include realistic treatment of 
gas and stellar feedback, which is important for the survival 
of a bar. Bars form relatively late: well after the last major 
merger (z ~ 1 — 2 for normal spiral such as our Milky Way), 
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when a collision of gas rich galaxies brings lots of gas with 
substantial angular momentum to the central disk galaxy. 
As the disk accretes the cold gas from the halo, it forms 
a new generation of stars and gets more massive. At some 
stage, the disk becomes massive enough to become unstable 
to bar formation. Once the stellar bar forms, it exists for the 
rest of the age of the Universe. 

The cosmological simulations are still in preliminary 
stages, and it is likely that many results will change as they 
become more accurate and the treatment of the stellar feed- 
back improves. However, existing cosmological simulations 
already show us the place and the role of traditional TV-body 
simulations of barred galaxies, which start with an unstable 
stellar disk. It was not clear whether and how this happens 
in the real Universe. Now the cosmological hydrodynamic 
simulations tell us that this is somewhat idealistic, but still 
a reasonable setup compatible with cosmological models. 

Simulations of bars play an important role for 
under standing th e phenomenon of barred galaxies 



(e.g., [M iller 



th e ph enomenon ot barred g a laxies 
19781; ISellwopdl Il980t lAthanassoulal 120031 ; 



IValenzuela fc Klypin I l2003h . Numerical models 



fully acc ount for many observed features of real barred 
galaxies jBureau fc Athanassoulall2005l ; iBureau et alj|2006l ; 
iBeaton et alj|2007h . So, it is very important to assess the ac- 
curacy of those simulations. Recent disagreements between 
resul ts of different research grou p s (IValenzuela fc Klypin 



20031 ; lO'Neill fc Dubinski I 120031 ; ISellwood fc Debattista 



2006) prompted us to undertake a careful testing of numer- 
ical effects and to compare results obtained with different 
codes. This type of code testing is routine in cosmologica l 
simulations (jFrenk et al.ll 19991 ; iHeitmann et al.ll2005l . l2007( l. 
but it never has been done before for bar dynamics. 
Testing and comparison of numerical codes is important for 
validating different numerical models. It was instrumental 
for development of precision cosmology. It is our goal to 
make such tests for TV-body models of barred galaxies. 

W e use four d i fferen t popular TV-body codes: 
' 19971), G a dget- 1 and Gadget-2 
Springel I 120051 ), and Pkdgrav 



|Sprine:el et al. 


|200l|; 


JWadslev et al. 


2004). 



ART is an Adaptive Mesh Refine- 
ment code that reaches high resolution by creating small 
cubic cells in areas of high density. Gadget and Pkdgrav, 
on the other hand, are TREE codes that compute forces 
directly for nearby particles and use a multipole expansion 
for distant ones. We use the codes to run a series of 
simulations using the same initial conditions for all codes. 

We also address another issue: the effects of disk thick- 
ness on the structure and evolution of bars. Only recently 
have the simulations started to have enough mass and force 
resolution to resolve the vertical height of stellar disks. We 
use different codes to show that the disk height plays an 
important and somewhat unexpected role. 

One of the contentious issues in the simulations of 
bar dynamics is the angular speed and the structure 
of bars in massive dark matter halos. The amount and 
the rate at which bars slow down is still under debate. 
iDebattista fc Sellwood ] (| 19981. l2000h find in their massive 
halo models, i.e., those for which the contributions of the 
disk and the halo to the mass in the central region are 
comparable, that the bar loses about 40% of its initial an- 
gular momentum, L z , in ~ 10 Gyr. However, in simula- 
tions with much better force resolution and a more real- 



istic cosmological halo setup, IValenzuela fc Klypin I (|2003T ) 
and lColm et al. I (120061) find a decrease in L z of o nly 4-8% 
in ~ 6 Gvr. IDebattista fc Sellwood I (|l998l . l2000h also find 
that bars do not significantly slow down for lower density 
haloS; 

IValenzuela fc Klypin I l|2003h presented bar models in 
which stellar disks were embedded in a CDM Milky- Way- 
type halos with realistic halo concentrations c ~ 15, where 
c is the ratio of the virial radius to the characteristic ra- 
dius of the dark matter halo. These simulations were run 
with the ART code with high spatial resolution of 20-40 
pc. The bars in the models were ro tating fast for billions 
of years. IValenzuela fc Klypin I (|2003h argued that slow bars 
in previous simulations were an artifact of low resolution. 
ISellwood fc Debattistal il2006h used initial conditions of one 
of the models of IValenzuela fc Klypin"! (120031 ) and run a se- 
ries of simulations using their hybrid, polar-grid code. They 
found, in most cases, a different evolution than that reported 
in Valenzuela & Klypin. In particular, contrary to Valen- 
zuela & Klypin's results, they did not find that the bar pat- 
tern speed is almost constant for a long period of time. They 
attribute t he differences to the ART refin ement scheme. 

While IValenzuela fc Klypin I (|2003h mention numer- 
ical effects (lack of force and mass resolution) as the 
cause for excessive slowing down of bars in earlier simula- 
tions, there was anot h er eff ect, which was not noticed by 
IValenzuela fc Klypin I (2003): Their disk was rather thin, 
with a scale-height h z of only 0.07 of the disk scale-length 
Rd. This should be compared with h z ~ 0.2 Rd used in 
most studies of stellar bars (e.g., Athanassoula fc Misiriotisl 
|2002| ; I Athanassoula 1120031; iMartinez-Valpuesta et al. 1120061 ). 
Models of IDebattista fc Sellwood I i|l99S!) have h, = 0.1R d , 
but the resolution of thei r simulations was grossly i nsuffi- 
cient to resolve this scale. lO'Neill fc Dubinski I (|2003l ) used 
h z = 0.1-Rd for a model, which had very little dark matter 
in the central disk region: Md m /Mdi s k ~ 1/4 inside radius 
R = 3R d- Dependence of bar speed on d isk thickness was 
noted bv lMisiriotis fc Athanassoula! (|2000l ): thicker disks re- 
sult in slower bars. 

The remainder of this paper is organized as follows. In 
section [2] we give a detailed review of available data on disk 
scale heights and present a simple analytical model for the 
relation of the disk scale length with the disk scale height. 
We describe our numerical models in Section [3] In Section [4] 
we give a brief description of the codes and present analysis 
of numerical effects. Main results are presented in Section [S] 
We summarize our results in Section [6] 



2 DISK HEIGHTS 

Disk scale height appears to play a very important role in 
the development of barred galaxies. Thus, it is important to 
know the range of disk scale heights in real spiral galaxies. 
The most accurate measurement of the disk height comes 
from the Milky Way galaxy. Edge-on galaxies provide an- 
other opportunity to measure disk heights, but those mea- 
surements are much less accurate because of dust absorp- 
tion close to the disk plane. There is not much disagree- 
ment between different studies regarding the disk thickness 
of the Milky Way: the ex ponential scale height of the stellar 
thin disk is h z « 300 pc (|Gilmore fc Reidlll983l ; lOiha et all 
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1 19991 ; | Juric et all 120081 ) . The thick disk component has a 
scale height of ~ 1 kpc, but it has a small fraction of mass. 
Neutral hydrogen and molecular gas have scales 200 pc and 
50 pc. So, we can estimate the scale height of the mass distri- 
bution as 250-300 p c. Using the exponential di s k scale length 
Rd y 3 kpc (e.g, iDehnen fc Binnevl Il998al ; iKlypin et al.l 
2002), we get the ratio of the scale height to the scale length 
hz/Rd «0.1. 

Observations of edge-on galaxies can be used 
to esti mate the disk heigh ts for ot her gala xies 

2002: 



van der Kruit fc Searld 
iBizvaev fc Mitronoval 120021: 



Yoachim fc Dalcanton 20061 ) 



1981: 
iKregel et al, 
Kreeel et all (|2005h 



Kreeel et al.l 
2005; 



gave 

ratios of exponential heights to exponential lengths for 34 
edge-on galaxies measured in 7 -band and f ound that the 
median ratio is h z /Rd ~ 0.12. ISeth et all |2005l ) present 
fits for brightness profiles in if-band (2MASS images) for 
two Milky Way- type galaxies ATGC891 (V m » = 214 km/s) 
and iVGC4565 (V ma x = 227 km/s). They find the half-light 
height to exponential disk scale ratios zx/ijRd = 0.072 
and 0.085 for t he two galaxies correspondin gly. For the five 
galaxies in the lYoachim fc Dalcantonl l[2006j) sample, which 
had circular velocities in the range 150 — 200 km/s, the 
average ratio was Zi/ 2 /Rd ~ 0.1. 

Estimates of disk heights in edge-on galaxies suf- 
fer f rom substantial absorption close to the plane of the 
disk l|Xilouris et al.lll999l ; lYoachim fc Dalcanton|[2 006). This 
makes scale heights of the thin disk difficult to measure di- 
rectly and causes the results to be dominated by flux coming 
from high galactic latitudes, where the thick disk is domi- 
nant. In turn, this leads to a substantial overestimation (by 
a factor of 2-3) of the disk heights even in red bands, if one 
interprets those as estimate s of the thin disk component 
|Yoachim fc Dalcantonll2006l ). 

We can use stellar-dynamical arguments to estimate the 
scale-heights. The idea is to use the ratio of the vertical 
velocity dispersion a z to the radial velocity dispersion or 
l|Kregel et al.|[2005l) . For old main sequence stars (B — V > 
0.5) in the solar nei ghborhood, the ratio is m easured to be 
Oz/or = 0.5 - 0.6 (jPehnen fc Binnevl Il998bl) . We find the 
same ratio for most of our dynamical model. The vertical 
velocity dispersion is related to the disk scale-height, and 
the radial velocity dispersion is related to the disk scale- 
length (among other parameters of the disk). 

Assuming an exponential stellar disk with the vertical 
density profile sech 2 (z/h z ), one gets: 

3.36G£(i?) 



ctr(R) = Q- 



k{R) 



a 2 (R) =TvGh z T,(R), 



(1) 



where Q is the Toomre stability parameter; £(i?) is the sur- 
face density, and k(R) is the epicycle frequency. For galaxies 
with fiat rotation curves and for radii R > Rd we can use 
k = V2V c i IC / R. The circular velocity V c i rc is defined by the 
mass distribution given by the sum of three components: 
disk, bulge, and dark matter. It is convenient to parameter- 
ize those relative to the total disk mass: M{R) /Mdisk = 
/di 8 k(-R) + M bulgc /M disk + M Am {R) /M disk , where / disk (.R) 
is the fraction of the disk mass inside radius R. Combin- 
ing these relations, we get the following expression for the 
height-to-length ratio: 



h z 
Rd 



/a z 3.36Q V 
\<tr 2tv ) 



X 3 exp(-X) 



/disk (X) + /bulge + /dmpf) ' 



(2) 



where X = R/Rd, /bulge = A/buige/M disk , and fdm(X) = 

Afdm(A)/Af d isk, 

We can now estimate the disk height at different radii. 
For example, we can get it at R — 3Rd assuming that the 
mass of dark matte r is about equal to th e disk mass /dm = 1 
|Klypin et alj|2002l ; IWidrow et ai1)2008l) and taking a small 



bulge / b ui g e = 0.2. For Q = 1.5 l|Widrow et al.ll2008l ) and 
taking o z /(tr = 0.5 we get h z = 0.11-Rd, which is consistent 
with the height of Milky Way disk. Eq. ([2} gives nearly the 
same height for Rd < R < 3Rd, if we scale the dark matter 
contribution in such a way that the rotation curve stays 
constant. 

To summarize, the disk scale-heights are relatively small 
for high surface brightness galaxies such as our Milky Way 
with h z = 0.1-Rd being a reasonably accurate estimate. Care 
should be taken not to over interpret results from edge-on 
galaxies. 



3 MODELS AND SIMULATIONS 
3.1 Initial conditions 

The setup of initial condit ions is described in detail in 
IValenzuela fc Klypin 1 12003). Here we briefly summarize the 
most important features. The system of a halo and a disk, 
with no initial bu lge or bar, is generated using the method of 
iHernquist I ljl993l ). In cylindrical coordinates the density of 
the stellar disk is approximated by the following expression: 



Pd(R,z) = 



M d 



4Trh z R 2 d 



sech' 



(ft, 



(3) 



where Md is the mass of the disk, Rd is the scale length, 
and h z is the disk scale height. The latter is assumed to be 
constant through the disk. The radial and vertical velocity 
dispersion are given by eqs. (JTJ. Our models keep Q fixed 
along the disk. The azimuthal velocity and its dispersion are 
found using the asymmetric drift and the epicycle approxi- 
mations. 

The mod els assume a NFW density profile 

|Navarro et al. I Il997l ) for the dark matter halo com- 
ponent, which is described by 



Pdm{t) = 



po 



M vir = 4-Kp r s 



x(l + x) 2 
m(l + c) 



1 + c 



r/r s , 

Rvir 



(4) 
(•>) 



where M v i r , Rvir, and c are the virial mass, the virial radius, 
and concentration of the halo, respectively. Given Mvir, the 
virial radius i s found once a cosmo l ogy is adoptecQ. Equa- 
tions (4-56) of iBinnev fc Tremaine I (1 19871 ) and the assump- 
tion of isotropy in the velocities allow us to determine the 
radial velocity dispersion as 



>r,DM — 



PDM 



GMM , 
Pdm 5 — dr, 



(6) 



where M(r) is the mass contained within radius r and G the 
gravitational constant. 



1 We adopt the flat cosmological model with a non-vanishing 
cosmological constant with Qq = 0.3 and h = 0.7. 
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Table 1. Initial Parameters of the models 



Code 


Name 


N disk 






Force resolution 


Time-step 


Disk scale height h z 






(10 s ) 


(io 6 ) 


(10 ) 


(PC) 


(ioV) 


(P c ) 


(!) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 




K series models: Nldisk = 


5 x 10 10 


M Q Utot = 1.43 x 1O 12 M R d 


= 3.86 kpc Q = 1.8 c= 10 


ART 


K al 


2.33 


2.7 


6.6 


14 


1.4 


200 


ART 


K a2 


2.33 


2.3 


6.6 


86 


1.9 


200 


ART 


K a3 


2.00 


2.2 


5.9 


170 


2.2 


714 


Gadget- 1 


Kgl 


1.00 


1.1 


2.9 


280 


8.6 


714 


Gadget-2 


Kg2 


2.33 


2.5 


6.7 


112 


2.6 


200 


Gadget-2 


K 9 3 


4.67 


5.0 


13.8 


112 


3.3 


200 


Gadget-2 


Kgi 


2.00 


2.2 


5.9 


140 


1.4 


200 


Gadget-2 


K 95 


1.00 


1.1 


2.9 


280 


29.2 


714 


Gadget-2 


K g6 


1.00 


1.1 


2.9 


280 


3.6 


714 


Pkdgrav 


K p i 


1.00 


1.1 


2.9 


136 


24.5 


714 


Pkdgrav 


K p2 


2.33 


2.5 


6.7 


136 


1.2 


200 




Model D 


: M- disk 


= 5 x 10 


1O M M 


tot = 1.43 x 10 12 M 


Q R d = 2.57 kpc Q = 1.3 c = 17 


Gadget-2 




2.33 


2.5 


6.7 


112 


2.6 


200 



3.2 Description of the models 

Selection of parameters of our models is motivated by a 
number of reasons. First, to simplify the comparison with 
previous res ults, we chose para meters, which are close to 
those used in ICoh'n et al. I (|2006l ) . Indeed, some of our mod- 
els have_exactly the sa me parameter as models Khb and D cs 
in ICoh'n et al. I (|2006T ). In this paper we preserve the first 
letter of the model name (K or D), but use subscripts to 
identify numerical code. Second, in order to test the effects 
of the disk height, we construct a new model by taking the 
Khb model and giving it a larger scale-height h z = 714 pc 
instead of h z — 200 pc. Third, our models are motivated by 
predictions of cosmological models. Thus, the models have 
extended dark matter halos with the NFW profile, realis- 
tic virial masses, and concentrations. The central regions of 
the models are dominated by the disk. The radius at which 
initially the dark matter mass is equal to the disk mass is 
equal to 9 kpc for models K and to 6 kpc for models D. 
Initial profiles of different co mponents are prese nted in the 
top two panels of Figure 1 in ICoh'n et al.l (|2006l ). 

The disk scale-height increases in the course of evolu- 
tion. We find that for thin disk models it more than doubles 
after 5 Gyrs of evolution, while for disk thick models the 
scale-height increases less: by a factor of 1.6. As a result, 
the scale-height to scale-length ratio of evolved thin disk 
models is close to the observed h z /Rd = 0.1. The ratio for 
evolved thick disk models is, on the other hand, twice the 
observed one. 

Table 1 presents the parameters of the models. The first 
and the second columns give the name of the code and the 
name of the model. The capital letter of model name rep- 
resents the model type (K or D). The first subscript in the 
model name indicates the code used to make the run: a for 
ART, g for Gadget, and p for Pkdgav. In columns (3), (4), 
and (5) we show the disk, the total (disk + dark matter), and 
the effective number of particles - the number of particles, 
which we would need, if we used equal-mass particles. The 
force resolution (6-th column) is twice the smallest cell size 
for ART code and the spline softening (2.8 times the effective 
Plummer softening) for tree codes. Note that the force reso- 



lution is the distance at which the force accurately matches 
the Newtonian force. The force continues to increase even 
below the resolution resulting in large changes in density. 
Column 7 shows the smallest time-step of simulations. All 
codes use variable time-steps. Details are given in the next 
section. 

The number of particles inside a sphere of radius equal 
to the force resolution is quite substantial. For a typical 
simulation such as K a 2 or K g 2 with the resolution ~ 100 pc, 
there are ~ 100 particles inside the resolution radius at the 
center of the system at the initial moment. For TREE codes, 
which keep the resolution constant, the number declines with 
distance, but it is still large in the plane of the disk. For 
example, it is ~ 10 — 15 at 8 kpc. The ART code maintains 
the nearly constant number of particles inside (increasing) 
radius of the force resolution at the level of 60 particles (see 
details in lColm et al~ll2006h . 

Particles with different masses are used in our simula- 
tions to increase the mass and force resolution in the central 
disk region. This is done by placing many small-mass par- 
ticles in the central disk-dominated region and by utilizing 
large-mass particles in the outer halo-dominated areas. We 
use four mass species. The first species represents disk and 
dark matter particles in the central halo region (the central 
~ 40 kpc region). Both the disk particles and the central 
dark matter particles have the same mass. More massive 
particles rarely enter the central ~ 10 kpc region. The mass 
species differ by a factor of two between one species and the 
next. 

The mass resolution - the mass of a disk particle or the 
mass of the smallest dark matter particle - is given by the 
ratio of the disk mass to the number of disk particles. It is 
in the range mi = (1 — 5) x 1O 5 M0. We present the time- 
step and the initial scale height in the last two columns. The 
typical number of time-steps for simulations is (2 — 4) x 10 5 . 
The physical parameters of the models, such as the mass of 
the disk or concentration of the dark matter halo, are shown 
in separate rows. 

In order to estimate the bar pattern speed Q p we first 
determine the orientation of the bar by iteratively applying 
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t(Gyr) 

Figure 1. Effects of the time-step on the evolution of bars. The 
time-step for each model is indicated in the right corner of the 
top panel. Evolution of the pattern speed of the bar (top panel) 
and the bar amplitude (lower panel) is shown for a thin disk K 
model with force resolution 560 pc. Models only start to deviate 
from each other at the beginning of the buckling phase (t f» 2.5 — 
3 Gyrs), which is marked by a slight drop in the amplitude of 
the bar. Once this stage is passed models with sufficiently small 
time-steps start to converge again except for the run with the 
largest time-step. We find that the time-step dt s=s 10 5 yrs is not 
sufficiently small for this model and produces incorrect results. 



the method of the inertia tensor in the plane of the disk. 
Q p is obtained subsequently by numerical differentiation: 
Q p = d<f>/dt, where 4> is the position angle of the bar. fn prac- 
tice, we use about ten consecutive snapshots for which the 
increasing function <f) is available and make a least squares 
fit. Then Q p is given by the slope of th e straight line. The bar 
ampli tude A-2 is computed similar to IValenzuela fc Klypin I 
(2003). For each logarithmically spaced cylindrical bin we 
find the amplitude of the second Fourier harmonic. Then the 
amplitude is smoothed over the radius and the maximum is 
taken as the bar amplitude A 2 . The bar length is defined 
as out-most radius at which the ratio of axes of isosurface 
density contours is 1.5. When finding radius of corotation, 
we use rotational velocity curve of disk particles Vrot(r) and 
find radius, at which it is equal to fl p r. 



4 NUMERICAL EFFECTS 
4.1 Codes 

Simulations were run w ith three TV-body code s: ART (Adap 
tive Refin ement Tree, Kravtsov et alj 19971). 



ment (AMR) TV-body code, which achieves high spatial res- 
olution by refining the base uniform grid in all high-density 
regions with an automated refinement algorithm. Gadget is 
a parallel TREE code. Here we only use the TV-body part of 
the code. Pkdgrav is another TREE code. These two TREE 
codes differ in the characteristics of the gravitational tree 
algorithm. For instan ce, Gadget employs th e Barnes & Hut 
TREE construction dBarnes fc Hut I I 1986ft while Pkdgrav 
uses a binary K-D TREE (iBentlev 1 1 1979ft . Gadget-2 only 
uses monopole moments in the multipole expansion while 
Pkdgrav advocates a hexadecapole as the optimum choice. 
(Gadget-1 can be configured to use octopole moments). The 
codes also differ in the cell-opening criterion. Here we use an 
opening angle criterion of 6 — 0.7 for Pkdgrav. In the case 
of Gadget-2 ru ns we use a tole rance parameter a = 0.005 
(see eq.(18) in ISpringel I (120051 )1 . The main advantage of 
Gadget-2 as compare d with Gadget- 1 is a more accurate 
time-stepping scheme (|Springel 112005ft . 

The codes use different variants of leapfrog scheme. An 
algorithm of integration of trajectories can be written as 
a sequence of operators, which advance particle positions 
(called drifts) and changes velocities (called kicks). For ex- 
ample, a simple constant-step leapfrog scheme is: 



v(t„+i/2) = v(t n _i/ 2 ) + g(t„)dt Kick, 
x(t n+1 ) = x(t n ) + v(t n+1 / 2 )dt Drift. 



(7) 
(8) 



Thus, the leapfrog integration is a sequence of 
K(dt)D(dt)K(dt)D(dt)... operators. We call this a KD 
scheme. Note that the order of operators is not important: 
the DK scheme is identical to the KD scheme. When the 
time-step changes with time, the order of operators makes 
a difference. It is also important to select the moment, 
at which t h e tim e-step should be changed. Following 
iQuinn et al l lll997h we use the operator S to indicate 
the moment when a new time-step is selected. We also 
need to specify the time when accelerations (in kicks) and 
velocities (in drifts) are estimated relative to the moment 
to which the velocities and coordinates are advanced. 
We attach the sign '+' ('-') to the name of operator 
to indicate that the operator uses information from the 
beginning (end) of the time-step. For example, K-(dt) 
is: v(t„) = v(tn-i) + g(t n -\)dt and the operator D+(dt) 
is x(t n ) = x(t n ~i) + v(t n )dt. Operators with subscript 
are time symmetric: they use information at the middle 
of the time-step. Using these operators we can write the 
algorithms of time-stepping in all our TV-body codes: 



SD-(dt/2)Ko(dt)D+(dt/2), Gadget - 1 

SK-(dt/2)D {dt)K+{dt/2), Gadget - 2 

K+(dt 1 /2)SK-(dt/2)D (dt), ART 

D-(dt/2)SK {dt)D+(dt/2). Quinn et al. 



(9) 
(10) 

(11) 
(12) 



Gadg et-2 dS pringcl ct al. 2001|; ISpringel 



I~[ |2 

grav l|Wadslev et al. 112004 ). ART is an adaptive mesh refine- 



Gadget-1 or 
2005ft , and Pkd- 



The Pkdgrav code has different integration schemes. 
T he scheme given in e qs. (|12jl uses the algorithm described 
bv lQuinn et al.l l|l997h . In our simulations with the Pkdgrav 
code, we use the scheme which is identical to the Gadget-2 
code both in the sense of the sequence of stepping and refine- 
ment conditions . This is, in fact, the way Pkdgrav is most 
commonly run dWadslev et al. Iliooi ). The ART code uses 
the time-step dti from previous moment to start the inte- 
gration. Then, having information on coordinates, it makes 
a decision on the value of new time-step. The Gadget-2 and 
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Figure 2. Errors of orbit integration. We integrate trajectory of a 
particle moving in a gravitational potential created by isothermal 
distribution of matter c/> = log(r), p tx r~ 2 . The particle has 
angular momentum 1/5 of the circular orbit and moves on an 
elongated orbit with r max /r m ; n m 10, which is not unusual for 
orbits in strong bars. The number of time-steps is 200 per radial 
period, which corresponds to the time-step of 10 5 yrs when scaled 
to a realistic galaxy model. The top panels show the error in 
energy conservation and the bottom ones are the errors in position 
angle. The left column is for a the leap-frog scheme implemented 
in Gadget-2 and Pkdgrav codes. The middle (right) column is for 
ART (Quinn et al.) code. The energy conservation is reasonably 
small, but errors in the orbital angle are unacceptably large. 



the ART schemes look different, but actually they are iden- 
tical, which can be seen when one writes the sequence of a 
few time-steps. The Gadget-1 and the Quinn et al. schemes 
look very similar, but they are quite different: the position of 
the S operator mak es the Quinn et al. scheme more accurate 
l|Quinn et al.lll997T ). 

Conditions for changing the time-step are different in 
different codes. In the ART code the time-step decreases by 
factor 2 when the number of particles exceeds some speci- 
fied level (typically 2-4 particles). A cell that exceeds this 
level is split into eight smaller cells resulting in the drop by 
2 3 times of the number of particles in a cell. This prescrip- 
tion gives scali ng of the time-step with the local density p 
as dt tx p -1 ^ 3 . IColfn et al. I (|2006T ) give more details of the 
procedure. The time-step in Q uinn et al. scheme scales as 
dt tx o~ 1 ^ 2 . [Zemp et all l|2007l ) also advocate a scheme with 
this scaling of the time-step. The Gadget and Pkdgrav code 
use a scaling with the gravitational acceleration dt tx g~ x ^ , 
which for p tx r~ 2 gives dt tx p _1,/4 . Among all the codes 
the Quinn et al. scheme uses the most aggressive prescrip- 
tion for changing the time-step and Gadget has the smallest 
change in dt. 

In practice, the time-step dt for the Gadget-2 code is 
defined by parameter 



Figure 3. Errors of orbit integration. The same as in Figure [2] 
but for 2000 time-steps per radial period. This corresponds to a 
time-step ~ 10 4 yrs in realistic simulations. The errors of integra- 
tion are very small for both the energy conservation and for the 
orbital angle. 



where |g| is the acceleration, and e is the effective Plum- 
mer softening. For example, we used rj — 0.04 for D g i and 
K g 2 models. It is T) = 0.11 for K g i model. In simulations 
where we increase the time-step by factor two, we double rj. 
For the Pkdgrav code the time-step is defined by a similar 
expression: if = |g|dt 2 /e s , where e s is the spline softening. 
We use rj = 0.025 for model K P 2 and r\ = 0.050 for model 



K 



pi- 



2 

7/ = 



\s\dt 2 
2e 



(13) 



In all codes the time-step changes discretely by factor 
two: dt — dto2~ m , where dto is the maximum time-step and 
m is an integer defined by local conditions (local density or 
local acceleration). 



4.2 Numerical Effects: resolution and 
time-stepping 

In order to obtain accurate results the trajectories of parti- 
cles should be integrated with a sufficient precision. There 
is no reliable method of estimating how small a time-step 
should be. While there are theoretical argumen ts what inte- 
grati o n schemes should (or should not) be used dQuinn et all 
Il997t iPreto fc Tremainei Il999l : ISpringel I |2005T> . only tests 
can tell how accurate results are. Simulations of bars have 
a special reason why the orbits should be accurate. Parti- 
cles, which make up the bar move typically on quite elon- 
gated trajectories periodically coming close to the center. 
When this happens, the acceleration changes substantially, 
and fast changing accelerations pose problems for numeri- 
cal integration. If accuracy of integration is not sufficient, a 
particle may erroneously change its direction of motion and 
start moving away from the bar. This artificial scattering on 
the center results in a smaller number of particles staying 
in the bar. Thus, it is important to have not only accurate 
particle energies, but also to have accurate phases of trajec- 
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Figure 5. The evolution of the bar pattern speed f2 p (top panel) and the bar amplitude A2 (lower panel) for K thick-disk models. All 
these models evolve similarly: it takes ~ 2 — 3 Gyrs to form a bar; the buckling phase happens at t ~ 5 Gyrs followed by a regime of a 
constant amplitude and nearly constant pattern speed. Through the course of evolution the pattern speed declines by a factor 2-3 and 
it is very low. 



tories. The later is more difficult than the former because 
different leap-frog schemes used in current N— body codes 
are known t o have problems with accurately tracking orbital 
phases (e.g.. ISpringel ||2005| ). 

The problem is more complicated because there is a real 
scattering on a dense central region: trajectories have a ten- 
dency to get deflected, when they come close to the center. 
The magnitude of this effect depends on the mass and size 
of the central concentration. Here the force and mass reso- 
lution play important role. If the resolution is not sufficient, 
the density in the central region will be lower resulting in 
less deflection of orbits. Thus, more particles will stay in the 
bar, which starts to trap even more particles. This leads to 
excessive growth of the bar. 

In order to estimate the effects of time integration of or- 
bits, we start by making simulations of the same model using 
different time-steps. We use the Gadget-2 code to make sim- 
ulations of the model K with low force resolution of 560 pc 
and with Ndisk = 10 5 . The initial disk height is h z = 200 pc. 



Figure \T\ shows results for simulations with the time-step 
changing by factor two from one simulation the other. Over- 
all, we clearly see convergence of the results, but it is not 
monotonic. Runs with different time-steps evolve very sim- 
ilarly until (2.5-3) Gyrs, when they start to diverge. The 
reason for the divergence is likely to be related with the fact 
that at around 2-3 Gyrs the system goes through the buck- 
ling instability. At this stage the errors in orbit integration 
produce large errors in the system configuration. For this 
particular model the time-step of 10 yrs was not sufficient: 
it resulted in a qualitatively wrong answer. 

It is difficult to predict what actually happens when 
the time-step is not small enough. In this case it produced 
a steep decline in the bar amplitude during the buckling 
stage. We had the same effect (weakened bar) in ART sim- 
ulations done with too large time-steps. At the same time, 
it is quite possible to have just an opposite effect: an ar- 
tificia lly stronger bar with l o wer p attern speed. For exam- 
ple, ISellwood fc Debattista I ([2006) run a model similar to 
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Figure 4. Dependence of bar properties on the force resolution. 
Evolution of the pattern speed of the bar (top panel) and the 
bar amplitude (lower panel) are shown for the high resolution 
model K g 2 with force resolution 112 pc)( full curves) and for a 
low resolution simulation with 560 pc resolution (dashed curves). 
The later simulation is the small time-step model presented in 
Figure[T]with the full curves. Increased force and mass resolution 
produce a more concentrated bulge, which weakened the bar. A 
shorter and weaker bar rotates faster and does not slow down 
much over many billions of years. 

K g 2 with a grossly insufficient time-step 1.5 x 10 5 yrs; 
they got a very long and a slow bar. For the same model 
IValenzuela fc Klypin I l|2003l ) used a ten times smaller step 
and got a significantly shorter and faster bar. 

The accuracy of orbit integration greatly depends on 
the distribution of density in the central region: the steeper 
is the density profile, the more difficult is the simulation. 
In our case the models K with a thick disk (h z = 714 pc; 
e.g., runs K g s and Ka.3) did not form a dense center. We 
find that for these models a relatively large step of dt ~ 
10 5 yrs is sufficient. Models K with thin disk (h z = 200 pc) 
produce a nearly flat circular velocity curve implying a steep 
profile, which can be roughly approximated as p oc r~ 2 in 
the central 2 kpc region. These models showed inconsistent 
results when large time-steps dt > 10 5 yrs were used. Only 
when we changed to much smaller steps {dt ~ 10 4 yrs), 
results became stable. 

The observed effect of time-stepping presented in Fig- 
ure[T]is somewhat unexpected. The most difficult and impor- 
tant part of the simulation is the motion of particles in the 
central region. Taking a typical particle velocity of 200 km/s 
we estimate that it would take 2 x 10 7 yrs and 200 time-steps 
(dt = 10 5 yrs) for a particle to cross the central 2 kpc re- 
gion. One would expect that 200 time-steps is enough to 
provide reasonably accurate results. Unfortunately, this is 



not the case as Figure \T\ shows. Indeed, our simple tests 
described below indicate that in spite of the fact that the 
energy is reasonably well preserved, the orbital angle is not 
- trajectories are scattered very substantially when a large 
time-step is used. This numerical scattering may result in 
incorrect properties of the bar. 

We investigate the situation with the accuracy of or- 
bit integration by studying the motion of particles in an 
idealized, but realistic case of a spherical logarithmic poten- 
tial <j)(r) — log(r). We implemented four time-integration 
schemes used in our simulations: a constant leapfrog scheme, 
and block-schemes with a variable time-step used in ART, 
Gadget-2/Pkdgrav, and Quinn et al. codes. 

Figure [2] shows the results of integration of an eccentric 
orbit with the apocenter/pericenter ratio of about 10:1. The 
gravitational potential is normalized in such a way that that 
the binding energy is equal to the time averaged kinetic en- 
ergy. In order to find the accuracy of the orbital angle, we 
run the orbit with a very small time-step: 10 5 steps per pe- 
riod. Then we find the position angle of a low-accuracy run 
and compare it with the position angle at the same moment 
in the small-step run. As compared with a constant-step 
run, all variable-step integration schemes give smaller errors 
in the energy conservation, but the errors in the position 
angle are too large. The constant step run is definitely bet- 
ter, but even in this case the errors (~ 10°) are still large. 
Decreasing the time-step by factor of ten gives much better 
results as shown by Figure In this case even the phase of 
the orbit is accurately simulated. The difference between the 
ART and the Gadget codes is due to the fact that the ART 
code takes more small steps in the central part of the orbit 
(the total number of steps is the same). When we apply the 
Gadget time-step changing algorithm to the ART code, we 
get exactly the same results. Quinn et al. takes even more 
refinements in the center, which improves the accuracy. Yet, 
the code gives more accurate results even when it runs with 
the same time-step as Gadget. Still, for a 200 orbits integra- 
tion - typical for both the bar simulations and for existing 
cosmological runs - the differences between the codes are 
not essential. The Quinn et al. scheme wins when we make 
much longer runs: for 20,000 orbits with 200 step/orbit the 
Quinn et al. code gave 10 times better accuracy than ART, 
and GADGET was another factor of ten worse than ART. 

Adequate force and mass resolution is another impor- 
tant numerical issue, which is even more difficult to handle. 
The "converged" (in the sense of time-stepping) solution 
presented in Figure [T] has a strong bar, which dramatically 
slows down: the pattern speed declines by a factor of two 
over 5 Gyr period. Results presented in Figure|3]tell another 
story: a low-resolution simulation can be very misleading. 
In this case we use the same thin disk model with twice 
as many particles and run it with five times better force 
resolution. The high mass and force resolution simulation 
is qualitatively different: the bar is weaker and its pattern 
speed hardly changes. 



5 RESULTS 

Table 2 gives different properties of the simulated models as 
measured after 5 Gyrs of evolution. In column 2 we present 
the fraction of the angular momentum lost by the stellar 
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material. The pattern speed £l p and the ratio of corotation 
radius to the bar radius are presented in columns 3 and 4. 
The bar length is given in column 5. The last three columns 
give parameters of a double-exponential approximation of 
the stellar surface density: E(r) = Eb u i gc exp(— r/i?bui g e) + 
Edisk exp(-r/i?disk)- 

The models are clearly split into two groups: those, 
which started with a thick disk (K a 3, K g i, K g s, and K p i) 
and the models, which started with a thin disk. For exam- 
ple, the ratio of the corotation radius to the bar radius is 
about 1.7 — 1.8 for thick disk models. For thin disk models 
the ratio is visibly smaller: 1.2 — 1.4. Thin-disk models have 
shorter bars and less massive bulges. The differences are es- 
pecially striking when we compare simulations done with 
the same code and with similar parameters. For example, 
models K a 2 and K a 3 have very similar numerical param- 
eters (time-step, resolution, and number of particles). Yet, 
their parameters at 5 Gyrs (and actually at any moment) 
are drastically different. For example, the pattern speed for 
model K a 2 is 2.5 times larger than for the model K a 3- 



5.1 Thick Disk Models 

We start with analysis of thick disk K models. It takes 
~ 2 — 3 Gyrs to form a bar; the buckling phase happens 
at t ~ 5 Gyrs followed by a regime of a constant amplitude 
and nearly constant pattern speed. Figure [5] shows the evo- 
lution of the bar pattern speed Q p and the amplitude of the 
bar Ai for the models. In Figure [6] w e present the total cir- 
cular velocity profile -J GM(< r) /r, the radial and vertical 
stellar velocity dispersions, and the disk surface density as 
a function of distance to the center of the galaxy, R. The 
comparison is made at 6.5 Gyrs for three models: K a 3, K g s, 
and K p i. 

Figure [7] shows disk particles seen along the minor axis 
of the bar for K models with the thick disk K g $ (left panels) 
and K a 3 (right panels) at four different epochs. The selected 
time moments represent different stages of bar evolution. At 
4 Gyrs the bar is the strongest, and it has not yet started the 
buckling stage: the disk is still thin. At 5 Gyrs the system 
goes through the buckling instability. At this moment the 
models look different. The differences die out as the buckling 
instability proceeds. Once the buckling stage is finished, the 
models again are very close. 

Overall, the models evolve very similarly. Some small 
differences are observed during the buckling phase. Yet, 
models tend to converge at the end of the evolution. The de- 
gree of agreement between different codes as demonstrated 
by Figure[6]is remarkable. Inside the central 5 kpc region the 
surface density of the disk deviates from model to model by 
only few percent. The vertical velocity dispersion over the 
whole disk deviates not more than 5 km/s. 

This agreement between models demonstrates that the 
results are code independent: if simulations are done with 
sufficiently small time-steps and with similar force and mass 
resolution, all codes produce nearly the same results. The 
results also show that there are no problems with any par- 
ticular code: all codes produce the same "answer". 
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Figure 6. Profiles of total circular velocity (top-left panel), stellar 
radial and vertical velocity dispersions (top- and bottom-right 
panel, respectively), and surface density (bottom-left panel) are 
shown in this four-panel plot for thick-disk K models at 6.5 Gyr. 
Here we have selected one model for each code. They agree with 
each other remarkably well. 



5.2 Thin Disk Models 

Figures [8] and [9] show the evolution of the bar amplitude 
and the bar pattern speed for thin disk K models. Again, 
the simulations behave very similarly, but this time they 
are very different from the thick disk models K. The thin 
disk models do not show any substantial evolution in the 
bar pattern speed and in the bar amplitude. Nevertheless, 
the agreement between the simulations is not as close as in 
the case of the thick disk models. It is important to note that 
there are no systematic differences between results obtained 
with different codes. The comparison of two runs presented 
in Figure [9] is especially striking. For example, the Gadget 
model K g 2 shows some small decline in the pattern speed 
and some oscillations in bar pattern speed and bar ampli- 
tude. The same evolution of Q p is demonstrated by the ART 
model K a i, which also has the oscillations with a slightly 
larger amplitude. The bar amplitudes for the two models 
are also very close over the whole simulated period of time. 
At the same time, models K a 2 (ART) and K g i (Gadget), 
which give very similar results, do not show any indication 
of slowing down of the bar (see Figure |H). 

We believe that this agreement between results pro- 
duced by different codes indicates that there are real physical 
differences between the models and that those differences are 
not of numerical origin. In other words, when started with 
the same physical initial conditions, the evolution produces 
different results. In the case of the thin disk K models the 
differences are on the level of ~ 10 percent after 6 Gyrs of 
evolution. 

What are the possible reasons for these variations be- 
tween the models? The iV-body problem is deterministic: 
initial conditions uniquely define the outcome of evolution. 
If one starts with identical initial conditions and uses a per- 
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Table 2. Parameters of the models after 5 Gyrs of evolution 



Name 


AL/L 


Pattern speed 


^cor /-^bar 


Bar length 


Disk scale 


Bulge scale 


Bulge/Total 






n b 




Rb 


length 


length fi bu i gc 






(%) 


(Gyr- 1 ) 




(kpc) 


(kpc) 


(kpc) 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 




Thick disk models 












K a3 


9.9 


10.0 


1.7 


10.8 


6.7 


1.07 


0.38 




8.3 


10.3 


1.7 


10.4 


7.0 


0.97 


0.32 


K g5 


10.3 


9.1 


1.7 


12.1 


6.6 


1.10 


0.36 


K 9 6 


8.2 


10.5 


1.5 


11.5 


6.5 


1.03 


0.35 


K pl 


8.5 


10.3 


1.8 


10.0 


6.4 


1.00 


0.33 




Thin disk models 












K al 


5.1 


21.0 


1.25 


6.7 


5.7 


0.73 


0.26 


K a2 


3.3 


25.4 


1.16 


6.5 


5.2 


0.54 


0.20 


K 9 2 


2.7 


19.8 


1.15 


7.5 


5.8 


0.63 


0.23 


K 93 


2.6 


21.7 


1.20 


6.7 


5.3 


0.64 


0.21 




2.1 


22.8 


1.22 


6.0 


5.3 


0.61 


0.21 


K p2 


3.1 


19.3 


1.23 


7.4 


5.5 


0.66 


0.22 


Dgl 


10.1 


19.5 


1.4 


7.6 


5.2 


0.43 


0.35 
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Figure 7. Distribution of disk particles seen along the bar minor 
axis at four different epochs: K g s model (left panels) and K a 3 
model (right panels). Only particles with \y\ < 1.0 kpc are shown. 
In order to have similar number of particles in left and right panels 
only half of the particles, randomly chosen, were used in model 
K a 3- At moment 5 Gyr a bar buckling is clearly seen. 



feet code, the results must be unique. Yet, this is too sim- 
plistic. An unstable system may have divergent evolution- 
ary tracks. A fluctuation slightly changes the system, but 
because of instabilities, the fluctuation grows, and the final 
answer is different as compared with the evolution of the 
system without the fluctuation. The barred stellar dynami- 
cal models have two stages when a system is unstable: the 
initial stage of formation of the bar and the buckling insta- 
bility. In addition, even in the quiet periods, the bar itself 
is an example of a potentially unstable system. For exam- 
ple, it could start trapping more particles resulting in even 



stronger bar, which traps even more particles. To large de- 
gree, this is exactly what happened with the bars in thick 
disk K models. The bars were expanding until all the disk 
became the bar. 

The source of fluctuations is an interesting issue. In the 
simulations the fluctuations are related with the initial ran- 
dom phases and amplitudes and with numerical inaccura- 
cies. Yet, we should not forget that our models are only ap- 
proximations to the reality. In real galaxies there is no short- 
age of fluctuations including satellite galaxies and molecular 
clouds to name the few. Whatever is the source of fluctua- 
tions, simulations of the same physical model evolve slightly 
different and produce slightly different results. This is ex- 
actly what the thin disk models did. 



5.3 Effects of the phase-space density: thin versus 
thick disks 

The only difference in the initial conditions between the thin 
disk and thick disk K models is the disk scale height. All the 
other parameters are the same. These seemingly minor vari- 
ations in initial conditions resulted in a remarkably differ- 
ent evolution and in large differences in the structure of the 
evolved models. In order to highlight the differences, Figure 
1101 shows the angular momentum loss AL / L of the stellar 
disk and the bar pattern speed for two models run with the 
Gadget code: one with a thin disk (Kgz) and another with a 
thick disk (-Kgs). The disk in the thick disk model loses four 
times more angular momentum and its bar rotates 2.5 times 
slower as compared with the thin disk model. The dark mat- 
ter is the same in both models, and it cannot be the reason 
of the difference. 

So, what are the possible reasons for such a large ef- 
fect of the disk height? The thickness o f the stellar disk af- 
fects the vertical waves and oscillations (|Merritt fc Sellwoodl 
Il994l ). This effect is definitely present and will have some 
impact on the evolution. At the same time, those vertical 
modes are likely to play a significant role during the buck- 
ling stage (|Merritt fc Sellwoodl H991 ). However, the differ- 
ences between the models develop too early and they are 
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Figure 8. The evolution of the bar pattern speed (top panel) and of the bar amplitude (lower panel) for thin-disk K models made with 
the Gadget-2, the ART, and the Pkdgrav codes. All the models show pattern speeds Q p , which do not show a systematic decline with 
time. There are differences between different runs, but they are relatively small: Cl p = (21.6±3)Gyr . Note substantial differences with 
respect to the thick-disk models presented in Figure [5] 



too large for the vertical modes to be the culprit. In the 
absence of a reliable theory of stellar bars, one can only 
speculate what is going on. One may think about two other 
effects, which can influence the systems: the Jeans mass and 
the phase-space density. The random velocities, which define 
the Jeans mass, tend to prevent collapse of perturbations 
on small scales. Thus, for the same random velocities, larger 
Jeans masses imply smaller densities. The phase-space den- 
sity acts in the same way. Thus, we expect that models with 
higher phase-space density (or small Jeans mass) will result 
in denser and more compact central region, which affects the 
growth of bars in a profound way. 

For a fixed surface density the vertical disk height h z 
changes the density in the disk p oc h~ x and the velocity 
dispersion o\ oc h z (see eq.((2}). Assuming for simplicity 
that the Jeans mass scales with the velocities and density 
as Mj oc oro^OzI yfp, we get Mj oc h z , which is a remark- 
ably strong effect considering that in our models h z varies 
by factor 3.5. The phase-space density changes even more: 



/ = p/(JR(J<t>Oz oc h z A ^ 2 . For our thin/thick models this gives 
variations by factor 7. 

The Jeans mass affects the evolution in a similar way 
as the Toomre stability parameter Q: larger Q results in 
later formation of the bar and in less prominent spiral arms. 
Comparison of the initial stages of bar formation in Figures[5] 
and[H]shows the effect: larger h z (thus, larger a z ) results in 
much delayed formation of the bar. 

The phase-space density is a complicated quantity, 
which in practice is used in the form of a c oarse-grained 
phase-space density. lAvila-Reese et all l|2005l ) present de- 
tailed results on the evolution of the phase-space density 
during the formation and evolution of barred disk models. 
As the system evolves, the coarse-grained phase-space den- 
sity / = P I a ro z decreases. The phase-space density / 
can be considered as a measure of the degree of compress- 
ibility of a gravitational system: for given rms velocities it 
defines the real-space density: p = fana ( f ) a z . The larger is 
the phase-space density, the larger is p. In turn, the den- 
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Table 3. Comparison of Parameters of the Milky Way and the model K g z 



Parameter K g $ 



Circular velocity (km/s) 220 

Surface disk density at Rq (M@/pc 2 ) 44.6 

Vertical rms velocity of stars at Rq (km/s) 14 

Radial rms velocity of stars at Rq (km/s) 38 

Pattern Speed Q p (km/s/kpc) 50 

Bar length (kpc) 3.3 

Total mass inside 60 kpc (1O 11 M ) 5.5 

Total mass inside 100 kpc (10 n M Q ) 7.3 



_ 40 ~ 




2 4 

Time (Gyrs) 



Figure 9. The evolution of the bar pattern speed (top panel) 
and of the bar amplitude (lower panel) for thin-disk models 
K g 2 (Gadget-2) and K a \ (ART). These two simulations were 
performed with different codes, yet they produce very similar re- 
sults. 



Milky Way Reference 

210-230 Binnev fc Merrifield (1998 . Sec. 10.6) 

48 ± 9 Kuiiken fc Gilmore (1991) 

15-20 Dehnen fc Binnev f 1998b) 

35-40 Dehnen fc Binnev (1998b) 

53 ± 3 Dehnen (1999) 

3.0-3.5 Freudenreich (1998) 



4 ±0.7 Xue et al. (2008) 




Time (Gyrs) 



Figure 10. Effects of the initial phase-space density: the evolu- 
tion of the pattern speed (top panel) and the disk angular momen- 
tum loss (bottom panel) for model K g s (full curves) and model 
K g r> (dashed curves). These two models differ by the initial value 
of the scale height of the disk, which results in substantially differ- 
ent phase-space density. The thin model goes through the buck- 
ling instability earlier and has a smaller angular momentum loss. 



sity in the central region of a bar is an important factor 
because according to our results it moderates the growth of 
the bar and can prevent it from growing excessively. Thus, 
the phase-space density in the central disk region is a funda- 
mental parameter, which significantly affects the evolution 
of barred galaxies. 

Indeed, the buildup of mass in the high phase-space 
density models happens in the models. Figure [TT1 shows the 
total (stellar plus dark matter) circular velocity (top-panel) 
and the stellar surface density (bottom-panel) profiles after 4 
Gyr of evolution. We see that the model K g 2 (thin disk) has 
a larger inner circular velocity and a larger surface density 
as compared with the low phase-space density model K g s. 
The rms velocities (right panels) are nearly the same for the 
two models. Thus, the K g i model has a substantially larger 
(by a factor of four) phase-space density. It started with a 
larger / and it ended with a larger /. 



6 DISCUSSION AND CONCLUSIONS 

We find that numerical effects - the mass and force resolu- 
tion, and the time integration of trajectories - can signifi- 
cantly alter results of simulations. The time-step of integra- 
tion must be small enough to allow accurate integration of 
expected trajectories. Accuracy of energy conservation can 
give a misleading impression that the simulation is adequate, 
while it actually makes substantial errors in position of or- 
bits. Our experiments with realistic orbits in models with 
flat rotation curves indicate that even the best available in- 
tegration schemes require ~ 2000 time-steps per orbit. This 
requirement is valid for variable-step schemes with more 
steps required for a constant-step schemes. Numerical tests 
with full-scale dynamical models confirm the condition. This 
condition results in very small time-steps of dt ~ 10 4 yrs. 
If one uses dimensionless units denned by scaling G = 1, 
Afdisk = 1) Rd = 1) then the dimensionless time step should 
be dt w 5 x 10~ 4 . This should be compared with typically 
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Figure 11. Comparison of profiles of models with different initial 
phase-space density. The circular velocity (top-left panel), the 
stellar radial and vertical velocity dispersions (top- and bottom- 
right panel, respectively), and the surface density (bottom-left 
panel), as a function of radius for models K g 2 and K g §, thin and 
thick disk, respectively. The comparison is made at 4 Gyr. The 
thin disk model (dashed curves) evolves to a more concentrated 
structure. The central concentration explains why the bar pattern 
speed for this model is higher than in the thick disk model. 




Figure 12. The iso-contours of the surface stellar density for 
representative models. The thin disk models K g -j, and D g \ de- 
velop realistic bars with the ratio of corotation radius to the bar 
length TZ = 1.2 (left panel; R COT = 4.5 kpc) and TZ = 1.4 (right 
panel; i? C or = 6.1 kpc). The thick disk model K g § (bottom panel) 
has a bar, that covers the whole disk with the corotation radius 
9.3 kpc. Distances in the plot are given in kpc units. All the mod- 
els are re-scaled to have the disk scale length 3 kpc - the same as 
for the Milky Way galaxy. 




10 



5 
X (kpc) 



10 



Figure 13. The iso-contours of the surface stellar density for 
model K g 3 (top panel) and the Milky Way galaxy (bottom 
panel). The K g 3 model was re-scaled to have the evolved disk 
scale length 3 kpc, which is close to the scale length of our Galaxy. 
The b ottom panel shows one of the models from Frcudcnrcich 
lll998l , fig. 14, right panel). The model represents a multi pa- 
rameter fit to the COBE DIRBE maps of the near-infrared light 
coming from central regions of our Galaxy. The small circle shows 
position of the Sun. The K„z model reproduces the length and 
the flattening of the Milky Way bar. 



used dt « 0.01 (e.g.. lAthanassoula fc Misiriotisl [20021 ) . This 
time-step would be insufficient for integration of models with 
central mass concentration presented i n this paper [fj A con - 
stant time-step dt — 5x 10 5 yrs used bv lWidrow et al.1 (|2008l ) 



2 Models studied by Athanassoula & Misiriotis (2002) are less 
concentrated and do not require a small time-step. 
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is too large for models of the Milky Way galaxy studied in 
their paper. 

The mass and related with it the force resolution also 
play important role. As Figure 2] shows, a low force resolu- 
tion produces a less dense central region, which results in 
a long and a very massive bar. The same model run with 
better force resolution produces a shorter, less massive, and 
faster rotating bar. 

Once the necessary numerical conditions are fulfilled, 
the codes produce practically the same results. We do not 
find any systematic deviations between the results obtained 
with TREE codes (Gadget and Pkdrgav) and with the 
Adaptive-Mesh-Refinement (ART) code. 

Disk height is an important parameter, which is often 
ignored in models of barred galaxies. In the models, which we 
consider in this paper, the disk height determines the global 
properties of the bars. Figure [12] shows surface density maps 
of the models with different initial disk thickness. Models 
with thin disks produce short bars with i?bar ~ -Rd, which 
rotate relatively fast: 1Z — 1.2 — 1.4 and which show very 
little decline of the pattern speed. Models with thick disk 
produce long and slow rotating bar. In order to facilitate 
the comparison with the our Galaxy, we re-scaled models 
to have the evolved disk scale length 2.65 kpc and to have 
the circular velocity at the solar distance 220 km/s. Any N- 
body system has two arbitrary scaling factors, which can be 
used to scale the system. 

Having scaled the models to fit the disk scale length, 
we can compare other parameters of the models. Because 
the simulations with thin disks produce reasonable models, 
we use one of the models {K g z) and compare it with the 
Milky Way. Table |3] gives a list of some parameters. Fig- 
ure [T3j£om£ares_the surface density maps of the Milky Way 
(|Freudenre ich 1998) and the K g 3 model. These comparisons 
show that the model fits the Milky Way reasonably well. 

We suggest that the disk height is only an indicator of 
a more fundamental property - the phase-space density in 
the central (R < Rd) region. In our models the phase-space 
density is uniquely related with the disk height. Initially 
in our models the disk height is low and the stability pa- 
rameter Q = 1.3 — 1.8 is constant across the disk. Thus, 
the phase-space density in the central region is high, and 
subsequent evolution brings that highly compressible stellar 
fluid close to the center, where it forms a nearly flat circu- 
lar velocity curve. The later, as we speculate, is responsible 
for arresting the growth of the bar. This relation between 
the disk height and the central phase-s pace density may not 
be tru e in general case. For ex ample, lO'Neill fc Dub inski 
(2003) and lWidrow et all {2008h consider models with large 
central Q and, thus, with a low central phase-space density. 
Indeed, in their models bars slow down substantially. 

The nearly constant pattern speed of bars in the thin 
disk models is somewhat puzzling. A bar is a very mas- 
sive non-axisymmetric object, which rotates inside a non- 
rotating dark matter halo. As such, one might expect that 
it should experi ence dynamical friction and s low down. This 
is why results o f lValenzuela fc Klypin I {2003 ) , which showed 
models with little slowing down of bars, were met with skep- 
ticism. Simxilati£ns_rjreserrte paper confirm the find- 
ings of lValenzuela fc Klypin I (2003) and show that they can- 
not be related with numerical problems. 

Orbital resonances may be responsible for the observed 



slow dynamical friction. The resonances in barred galax- 
ies have been extensivel y studied in recent y e ars using N- 
body simulations (e.g.. l Athanassoula I 120031: IColfn et al. I 
120061 ; ICeverino fc Klvpinll2007l ; I Weinberg fc Katzll2007l ). It 
is now well established that a large fraction of stellar par- 
ticles are in resonance with the bar. So me fraction of dark 
matter particles are also in resonance l|Coh'n et al. 1 120061 ; 
ICeverino fc Klypinl feoO?), and these are very important for 
the dynamical friction between the stellar bar and the dark 
matter. 

There are two types of (exact) resonances in 
an autonomous H amiltonian dynamical system: ellip- 
tic and hyperbolic (I Arnold fc Avezl Hi^S), Orbits, which 
are close to an elliptic resonance librate around the 
exact resonanc e and have the structure o f a sim- 
ple pendulum dLichtenberg fc Lieberrria 3 Il983l . Sec.2.4), 
{Murray fc Dermottl 1 19991 . Sec. 8). Hyperbolic resonances 
are points on intersections of separatrixes, which separate 
domains of elliptical resonances. Orbits close to a hyper- 
bolic resonance are unstable and tend to migrate away from 
the resonance, while orbits close to an elliptical resonance 
are stable. 

ICeverino fc Klypinl (|2007h study in detail the resonant 
orbits in simulations of barred galaxies. It appears that the 
orbits belong to elliptical resonances. These orbits track 
their resonance: the orbits do not evolve if the resonance 
does not move, and they follow the resonances if it gradu- 
ally migrates in the phase space. Thus, the elliptical reso- 
nances trap the orbits: if an orbit for whatever reason hap- 
pens to appear in the domain of the resonance, it will have 
a tendency to stay with the resonance. This phenomenon 
is thought to b e common in Solar System dynamics (e.g., 
lMalhotra|[l993l) . Resonance trapping explains why the sim- 
ulations show maxima in the distribution of orbital frequen- 
cies at the positi o ns of resonances. 

I Colin et al. I {2006) investigate another aspect of the 
resonant interaction between the stellar bar and the dark 
matter. They found that the dark matter particles, which 
are in resonance with the stellar bar, themselves form a bar, 
which rotates with the same angular speed and has a very 
small lag angle (~ 10°) as compared with the stellar bar. 
ICoh'n et al. 1 {2006) argue that the interaction between the 
dark matter and the stellar bars is the main mechanism for 
the dynamical friction between the disk and the dark mat- 
ter. The near alignment of the dark matter and stellar bars 
means that their interaction is minimized by the resonances. 

These results indicate that the resonant interaction be- 
tween the stellar bar and the dark matter is mostly due to 
elliptical resonances, and, thus, has a tendency to minimize 
the transfer of the angular momentum from the disk to the 
dark matter. Following orbits in such resonances emphasizes 
the need for conservative time-steps in iV-body simulations. 
Also subtle changes in the underlying global potential could 
change the relative number of orbits in these resonances 
leading to disparate results for the slowing of the bar. 
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